#!/usr/bin/env python3.9
# This generates 'analogue' SSTV from a (preferable very high resolution)
# image. The output is just in greyscale, because, well, it looks awesome.
# It's slow, and Python uses up a ridiculous amount of memory, but it's
# just accessible code for playing with.
# --
# Licenced according to the GNU Public Licence version 3.0
# Chris Baird,, May 2022 <vk2cjb@gmail.com>

import math, struct, sys
from PIL import Image

# scanlen is number of samples it needs = 48000 * 195.854ms, etc.
# The "Daytona SSTV Paper" scanline timings seem a bit off... QSSTV is happier
# with these (and probably other decoders using audio sampling at 48k).
# fudge is the 'true' duration of a scanline in audio samples, at least, as
# QSSTV sees it. It's a pain. QSSTV is probably faulty.
# And is QSSTV's hamdling of PD240 broken?

# Original specs mentioned in K0HEO/G4IJE's announcement
# Mode Pixel Line   Sync Porch Frame(secs) Resolution VIS VIS+P
# PD 90 532 703040  20000 2080 89.989120   320 x 256   99  99
# PD120 190 508480  20000 2080 126.103040  640 x 496   95  95
# PD160 382 804416  20000 2080 160.883200  512 x 400   98  226
# PD180 286 754240  20000 2080 187.051520  640 x 496   96  96
# PD240 382 1000000 20000 2080 248.000000  640 x 496   97  225

specs = [
    # name   VIS scanlen fudge     modesize
    ( "PD50", 93,  4394, 18632.85, (320,256)),
    ( "PD90", 99,  8173, 33748.00, (320,256)),
    ("PD120", 95,  5838, 24408.65, (640,496)),
    ("PD160", 98,  9389, 38614.60, (512,400)),
    ("PD180", 96,  8787, 36206.05, (640,496)),
    ("PD240", 97, 11736, 48003.25, (640,496)), # decodes poorly in QSSTV
    ("PD290", 94, 10984, 44992.50, (800,616))]
modelist = [x[0] for x in specs]
sample_rate, ph, nframes, mode, outfile = 48000, 0.0, 0, 4, ""
greymode, showident, leftident, redident = False, False, False, False
audio, domandel, infile = bytearray(), False, ""

###

while len(sys.argv) > 1:
    sys.argv.pop(0)
    s = sys.argv[0]
    if s == "-g": greymode = True
    elif s == "-i": showident = True
    elif s == "-r": redident = True
    elif s == "-l": leftident = True
    elif s == "-m":
        domandel = True
        greymode = True
        if not outfile: outfile = "mandelbrot"
    elif s == "-o":
        outfile = sys.argv[1]
        sys.argv.pop(0)
    elif s in modelist: mode = modelist.index(s)
    else:
        infile = s
        if not len(infile):
            print ("You needed instructions? oh.")
            sys.exit(1)

(name, VIS, scanlen, fudge, modesize) = specs[mode]
height = modesize[1]

#

if len(infile):
    im_orig = Image.open(infile)
    width, oheight = im_orig.width, im_orig.height
    if im_orig.mode == "I": depth = 65536.0
    else: depth = 256.0

    ycbcr = im_orig.convert('YCbCr')
    y, cb, cr = ycbcr.split()

    imY = y.resize((width,height)).convert("F")
    imCb = cb.resize((width,height//2)).convert("F")
    imCr = cr.resize((width,height//2)).convert("F")
else:
    width, oheight, depth = modesize[0] * 32, modesize[1], 65536.0

wantedwidth = (modesize[0] * oheight) / modesize[1]
scanlinescale = (modesize[0] * oheight) / (modesize[1] * scanlen)
leftedge = (width - wantedwidth) / 2

###

def sampleout(freq):
    global ph,nframes,audio
    s = int(math.sin(ph / sample_rate) * 3276.8) # constant can vary
    if s < 0: s += 65536
    audio += bytes([s & 255,s >> 8])
    ph += math.tau * freq
    nframes += 1


def tonegen(freq,samples):
    for _ in range(samples): sampleout(freq)

###

def vis(c):
    tonegen(1900,14400)  # leader - 1900 Hz, 300 ms/570 cycles/14400
    tonegen(1200,480)    # break - 1200 Hz,  10 ms/12 cycles/480
    tonegen(1900,14400)  # leader again
    tonegen(1200,1440)   # startbit - 1200 Hz, 30 ms/36 cycles

    p = 0
    for i in (1,2,4,8,16,32,64,128):
        if i == 128 and (p & 1): c |= 128
        if (c & i):
            p += 1
            tonegen(1100,1440)
        else:
            tonegen(1300,1440)
    tonegen(1200,1440)          # stop bit - 30 ms

###

id = (0,0,0,0,1,0,1,0,0,1,0,1,0,0,1,1,1,0,0,1,1,1,0,0,0,0,1,0,0,1,1,1,0,0,0,0,
      0,0,0,0,1,0,1,0,0,1,0,1,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,1,0,1,0,0,0,0,
      0,0,0,0,1,0,1,0,0,1,1,0,0,0,1,1,1,0,0,1,0,0,0,0,0,0,1,0,0,1,1,1,0,0,0,0,
      0,0,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,0,1,0,0,0,0,
      0,0,0,0,0,1,0,0,0,1,0,1,0,0,1,1,1,0,0,1,1,1,0,0,1,1,1,0,0,1,1,1,0,0,0,0)

def doidentline(y,side):
    yt = 36 * (int((60*y / modesize[1])) % 10)
    for x in range(scanlen):
        xx = int((72*x / scanlen)) % 36
        i = xx + yt
        freq = 1900
        if i < len(id) and id[i]: freq = 2300
        if not side and x < (scanlen/2): freq = 1900
        if side and x > (scanlen/2): freq = 1900
        sampleout(freq)

###

def mandel(a,b,n):
    x, y = 0,0
    for i in range(n):
        xx = x*x - y*y + a;
        y = 2*x*y + b;
        x = xx;
        if x*x + y*y > 4: break
    return i


def domandelline(y):
    xscale = modesize[0] / (modesize[1] * scanlen)
    for x in range(scanlen):
        # the classic mandelbrot pose
        #x1 = 0.8 - (2.5 * x * xscale)
        #y1 = -1.25 + (2.5 * y / modesize[1])
        #d = 2300 - 145.0*math.log(1+mandel(x1, y1, 300))

        # my favourite..
        x1 = -0.55 + (0.000011 * x * xscale)
        y1 = -0.55 + (0.000011 * y / modesize[1])

        m = mandel(x1, y1, 1500)
        #d = 2300 - 150.0*math.log(1+m)
        #d = 2300 - 150*math.log(1+m)
        d = 2300 - 114.3 * (m & 7)
        #d = 2300 - 800 * (m & 1)

        if d < 1500: d = 1500
        sampleout(d)


###

def getvalue(image,sx,y):
    x = leftedge + (sx * scanlinescale)
    if x < 0 or x >= width: r = 0
    else: r = image.getpixel((x,y)) / depth
    return r


def doline(image,y):
    for x in range(scanlen):
        c = getvalue(image,x,y)
        freq = 1500 + (800.0 * c)
        sampleout(freq)

### Main

tonegen(0,24000)
vis(VIS)

identstart, f, t = height * 0.80, 0, 0

for y in range(0,height,2):
    print (y, end='\r')
    f = int(y*fudge*0.5 - t)
    onf = nframes
    tonegen(1200,960+f)
    tonegen(1500,96)

    if domandel: domandelline(y)
    else: doline(imY, y)

    if greymode:
        if showident:
            if redident:
                if y > identstart: doidentline(y,leftident)
                else: tonegen(1900,scanlen)
                tonegen(1900,scanlen)
            else:
                tonegen(1900,scanlen)
                if y > identstart: doidentline(y,leftident)
                else: tonegen(1900,scanlen)
        else:
            tonegen(1900,scanlen)
            tonegen(1900,scanlen)
    else:
        doline(imCr, y/2)
        doline(imCb, y/2)

    if domandel: domandelline(y+1)
    else: doline(imY, y+1)
    #print (y, y*fudge*0.5, t, nframes - onf, f)
    t += nframes - onf

tonegen(0,24000)

###

if not outfile: outfile = infile + "-" + name + ".wav"

wav = open(outfile,"wb")

s = struct.Struct(b"<ccccLccccccccLHHLLHHccccL")
wav.write(s.pack(b'R',b'I',b'F',b'F', # riff_header
                 2*len(audio)+44,     # wav_size
                 b'W',b'A',b'V',b'E', # wave_header
                 b'f',b'm',b't',b' ', # fmt_header
                 16,                  # fmt_chunk_size
                 1,                   # audio_format = PCM
                 1,                   # num_channels
                 sample_rate,
                 sample_rate*1*2,     # byte_rate
                 2,                   # sample_alignment
                 16,                  # bit_depth
                 b'd',b'a',b't',b'a',
                 2*len(audio)))       # data_bytes
wav.write(audio)                      # thanks to bytearray
wav.close()

###
